!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

recursive subroutine block_boundary (current)

   use parameters
   use block_define
   use interfaces

   implicit none

   type(block), pointer, intent(inout) :: current

   integer(4) :: i, i_block, level, i_b, i_c, i_p, sig
   real(8) :: lmt, denominator
   real(8), dimension(ng + 1) :: ro_bnd, mx_bnd, my_bnd, mz_bnd, bx_bnd, by_bnd, bz_bnd, en_bnd
   type(block), pointer :: child, parent, neigh_l, neigh_r
!----------------------------------------------------------------------------------------------------------------------------------------------|

   level = current%lv
   parent => current%pp

!----------------------------------------------------------------------------------------------------------------------------------------------|
!      left boundary
!----------------------------------------------------------------------------------------------------------------------------------------------|
   if (current%neigh_lv(1) .ge. 1) then
      neigh_l => current%neigh(1, level)%pblock
      if (associated (neigh_l)) then
         do i = 1, ng
            i_b = current%nx - ng * 2 + i
            current%ro(i) = neigh_l%ro(i_b)
            current%mx(i) = neigh_l%mx(i_b)
            current%my(i) = neigh_l%my(i_b)
            current%mz(i) = neigh_l%mz(i_b)
            current%bx(i) = neigh_l%bx(i_b)
            current%by(i) = neigh_l%by(i_b)
            current%bz(i) = neigh_l%bz(i_b)
            current%en(i) = neigh_l%en(i_b)
         enddo
      else
         neigh_l => current%neigh(1, level - 1)%pblock
         if (associated (neigh_l)) then
            do i = 1, ng
               i_b = neigh_l%nx - ng * 2 + i
               ro_bnd(i) = neigh_l%ro(i_b)
               mx_bnd(i) = neigh_l%mx(i_b)
               my_bnd(i) = neigh_l%my(i_b)
               mz_bnd(i) = neigh_l%mz(i_b)
               bx_bnd(i) = neigh_l%bx(i_b)
               by_bnd(i) = neigh_l%by(i_b)
               bz_bnd(i) = neigh_l%bz(i_b)
               en_bnd(i) = neigh_l%en(i_b)
            enddo
            ro_bnd(ng + 1) = parent%ro(ng + 1)
            mx_bnd(ng + 1) = parent%mx(ng + 1)
            my_bnd(ng + 1) = parent%my(ng + 1)
            mz_bnd(ng + 1) = parent%mz(ng + 1)
            bx_bnd(ng + 1) = parent%bx(ng + 1)
            by_bnd(ng + 1) = parent%by(ng + 1)
            bz_bnd(ng + 1) = parent%bz(ng + 1)
            en_bnd(ng + 1) = parent%en(ng + 1)
            do i = 1, ng
               i_p = (ng + i + 1) / 2
               sig = (-1) ** (i + ng)

               denominator = ro_bnd(i_p + 1) - ro_bnd(i_p - 1)
               if (denominator .eq. 0) denominator = eps
               lmt = max((ro_bnd(i_p + 1) - ro_bnd(i_p)) * (ro_bnd(i_p) - ro_bnd(i_p - 1)), 0.0d0) / denominator
               current%ro(i) = ro_bnd(i_p) + sig * lmt * 0.5d0

               denominator = mx_bnd(i_p + 1) - mx_bnd(i_p - 1)
               if (denominator .eq. 0) denominator = eps
               lmt = max((mx_bnd(i_p + 1) - mx_bnd(i_p)) * (mx_bnd(i_p) - mx_bnd(i_p - 1)), 0.0d0) / denominator
               current%mx(i) = mx_bnd(i_p) + sig * lmt * 0.5d0

               denominator = my_bnd(i_p + 1) - my_bnd(i_p - 1)
               if (denominator .eq. 0) denominator = eps
               lmt = max((my_bnd(i_p + 1) - my_bnd(i_p)) * (my_bnd(i_p) - my_bnd(i_p - 1)), 0.0d0) / denominator
               current%my(i) = my_bnd(i_p) + sig * lmt * 0.5d0

               denominator = mz_bnd(i_p + 1) - mz_bnd(i_p - 1)
               if (denominator .eq. 0) denominator = eps
               lmt = max((mz_bnd(i_p + 1) - mz_bnd(i_p)) * (mz_bnd(i_p) - mz_bnd(i_p - 1)), 0.0d0) / denominator
               current%mz(i) = mz_bnd(i_p) + sig * lmt * 0.5d0

               denominator = bx_bnd(i_p + 1) - bx_bnd(i_p - 1)
               if (denominator .eq. 0) denominator = eps
               lmt = max((bx_bnd(i_p + 1) - bx_bnd(i_p)) * (bx_bnd(i_p) - bx_bnd(i_p - 1)), 0.0d0) / denominator
               current%bx(i) = bx_bnd(i_p) + sig * lmt * 0.5d0

               denominator = by_bnd(i_p + 1) - by_bnd(i_p - 1)
               if (denominator .eq. 0) denominator = eps
               lmt = max((by_bnd(i_p + 1) - by_bnd(i_p)) * (by_bnd(i_p) - by_bnd(i_p - 1)), 0.0d0) / denominator
               current%by(i) = by_bnd(i_p) + sig * lmt * 0.5d0

               denominator = bz_bnd(i_p + 1) - bz_bnd(i_p - 1)
               if (denominator .eq. 0) denominator = eps
               lmt = max((bz_bnd(i_p + 1) - bz_bnd(i_p)) * (bz_bnd(i_p) - bz_bnd(i_p - 1)), 0.0d0) / denominator
               current%bz(i) = bz_bnd(i_p) + sig * lmt * 0.5d0

               denominator = en_bnd(i_p + 1) - en_bnd(i_p - 1)
               if (denominator .eq. 0) denominator = eps
               lmt = max((en_bnd(i_p + 1) - en_bnd(i_p)) * (en_bnd(i_p) - en_bnd(i_p - 1)), 0.0d0) / denominator
               current%en(i) = en_bnd(i_p) + sig * lmt * 0.5d0
            enddo
            neigh_l%ro(neigh_l%nx - ng) = neigh_l%ro(neigh_l%nx - ng) - min_dt / neigh_l%dx * (current%flux_bnd(1, 1) - neigh_l%flux_bnd(2, 1))
            neigh_l%mx(neigh_l%nx - ng) = neigh_l%mx(neigh_l%nx - ng) - min_dt / neigh_l%dx * (current%flux_bnd(1, 2) - neigh_l%flux_bnd(2, 2))
            neigh_l%my(neigh_l%nx - ng) = neigh_l%my(neigh_l%nx - ng) - min_dt / neigh_l%dx * (current%flux_bnd(1, 3) - neigh_l%flux_bnd(2, 3))
            neigh_l%mz(neigh_l%nx - ng) = neigh_l%mz(neigh_l%nx - ng) - min_dt / neigh_l%dx * (current%flux_bnd(1, 4) - neigh_l%flux_bnd(2, 4))
            neigh_l%bx(neigh_l%nx - ng) = neigh_l%bx(neigh_l%nx - ng) - min_dt / neigh_l%dx * (current%flux_bnd(1, 5) - neigh_l%flux_bnd(2, 5))
            neigh_l%by(neigh_l%nx - ng) = neigh_l%by(neigh_l%nx - ng) - min_dt / neigh_l%dx * (current%flux_bnd(1, 6) - neigh_l%flux_bnd(2, 6))
            neigh_l%bz(neigh_l%nx - ng) = neigh_l%bz(neigh_l%nx - ng) - min_dt / neigh_l%dx * (current%flux_bnd(1, 7) - neigh_l%flux_bnd(2, 7))
            neigh_l%en(neigh_l%nx - ng) = neigh_l%en(neigh_l%nx - ng) - min_dt / neigh_l%dx * (current%flux_bnd(1, 8) - neigh_l%flux_bnd(2, 8))
         else
            write (*, *) 'Error from subroutine block_boundary:'
            write (*, *) ' Current level has no left neighbors and it is not a boundary block.'
         endif
      endif
   else
      call bnd_free_x (current%ro, ng, current%nx, 'l')
      call bnd_free_x (current%mx, ng, current%nx, 'l')
      call bnd_free_x (current%my, ng, current%nx, 'l')
      call bnd_free_x (current%mz, ng, current%nx, 'l')
      call bnd_free_x (current%bx, ng, current%nx, 'l')
      call bnd_free_x (current%by, ng, current%nx, 'l')
      call bnd_free_x (current%bz, ng, current%nx, 'l')
      call bnd_free_x (current%en, ng, current%nx, 'l')
   endif

!----------------------------------------------------------------------------------------------------------------------------------------------|
!      right boundary
!----------------------------------------------------------------------------------------------------------------------------------------------|
   if (current%neigh_lv(2) .ge. 1) then
      neigh_r => current%neigh(2, level)%pblock
      if (associated (neigh_r)) then
         do i = 1, ng
            i_b = current%nx - ng + i
            current%ro(i_b) = neigh_r%ro(i + ng)
            current%mx(i_b) = neigh_r%mx(i + ng)
            current%my(i_b) = neigh_r%my(i + ng)
            current%mz(i_b) = neigh_r%mz(i + ng)
            current%bx(i_b) = neigh_r%bx(i + ng)
            current%by(i_b) = neigh_r%by(i + ng)
            current%bz(i_b) = neigh_r%bz(i + ng)
            current%en(i_b) = neigh_r%en(i + ng)
         enddo
      else
         neigh_r => current%neigh(2, level - 1)%pblock
         if (associated (neigh_r)) then
            do i = 2, ng + 1
               i_b = i + ng - 1
               ro_bnd(i) = neigh_r%ro(i_b)
               mx_bnd(i) = neigh_r%mx(i_b)
               my_bnd(i) = neigh_r%my(i_b)
               mz_bnd(i) = neigh_r%mz(i_b)
               bx_bnd(i) = neigh_r%bx(i_b)
               by_bnd(i) = neigh_r%by(i_b)
               bz_bnd(i) = neigh_r%bz(i_b)
               en_bnd(i) = neigh_r%en(i_b)
            enddo
            ro_bnd(1) = parent%ro(parent%nx - ng)
            mx_bnd(1) = parent%mx(parent%nx - ng)
            my_bnd(1) = parent%my(parent%nx - ng)
            mz_bnd(1) = parent%mz(parent%nx - ng)
            bx_bnd(1) = parent%bx(parent%nx - ng)
            by_bnd(1) = parent%by(parent%nx - ng)
            bz_bnd(1) = parent%bz(parent%nx - ng)
            en_bnd(1) = parent%en(parent%nx - ng)
            do i = 1, ng
               i_c = i + current%nx - ng
               i_p = (i + 1) / 2 + 1
               sig = (-1) ** i

               denominator = ro_bnd(i_p + 1) - ro_bnd(i_p - 1)
               if (denominator .eq. 0) denominator = eps
               lmt = max((ro_bnd(i_p + 1) - ro_bnd(i_p)) * (ro_bnd(i) - ro_bnd(i_p - 1)), 0.0d0) / denominator
               current%ro(i_c) = ro_bnd(i_p) + sig * lmt * 0.5d0

               denominator = mx_bnd(i_p + 1) - mx_bnd(i_p - 1)
               if (denominator .eq. 0) denominator = eps
               lmt = max((mx_bnd(i_p + 1) - mx_bnd(i_p)) * (mx_bnd(i) - mx_bnd(i_p - 1)), 0.0d0) / denominator
               current%mx(i_c) = mx_bnd(i_p) + sig * lmt * 0.5d0

               denominator = my_bnd(i_p + 1) - my_bnd(i_p - 1)
               if (denominator .eq. 0) denominator = eps
               lmt = max((my_bnd(i_p + 1) - my_bnd(i_p)) * (my_bnd(i) - my_bnd(i_p - 1)), 0.0d0) / denominator
               current%my(i_c) = my_bnd(i_p) + sig * lmt * 0.5d0

               denominator = mz_bnd(i_p + 1) - mz_bnd(i_p - 1)
               if (denominator .eq. 0) denominator = eps
               lmt = max((mz_bnd(i_p + 1) - mz_bnd(i_p)) * (mz_bnd(i) - mz_bnd(i_p - 1)), 0.0d0) / denominator
               current%mz(i_c) = mz_bnd(i_p) + sig * lmt * 0.5d0

               denominator = bx_bnd(i_p + 1) - bx_bnd(i_p - 1)
               if (denominator .eq. 0) denominator = eps
               lmt = max((bx_bnd(i_p + 1) - bx_bnd(i_p)) * (bx_bnd(i) - bx_bnd(i_p - 1)), 0.0d0) / denominator
               current%bx(i_c) = bx_bnd(i_p) + sig * lmt * 0.5d0

               denominator = by_bnd(i_p + 1) - by_bnd(i_p - 1)
               if (denominator .eq. 0) denominator = eps
               lmt = max((by_bnd(i_p + 1) - by_bnd(i_p)) * (by_bnd(i) - by_bnd(i_p - 1)), 0.0d0) / denominator
               current%by(i_c) = by_bnd(i_p) + sig * lmt * 0.5d0

               denominator = bz_bnd(i_p + 1) - bz_bnd(i_p - 1)
               if (denominator .eq. 0) denominator = eps
               lmt = max((bz_bnd(i_p + 1) - bz_bnd(i_p)) * (bz_bnd(i) - bz_bnd(i_p - 1)), 0.0d0) / denominator
               current%bz(i_c) = bz_bnd(i_p) + sig * lmt * 0.5d0

               denominator = en_bnd(i_p + 1) - en_bnd(i_p - 1)
               if (denominator .eq. 0) denominator = eps
               lmt = max((en_bnd(i_p + 1) - en_bnd(i_p)) * (en_bnd(i) - en_bnd(i_p - 1)), 0.0d0) / denominator
               current%en(i_c) = en_bnd(i_p) + sig * lmt * 0.5d0
            enddo
            neigh_r%ro(ng + 1) = neigh_r%ro(ng + 1) + min_dt / neigh_r%dx * (current%flux_bnd(2, 1) - neigh_r%flux_bnd(1, 1))
            neigh_r%mx(ng + 1) = neigh_r%mx(ng + 1) + min_dt / neigh_r%dx * (current%flux_bnd(2, 2) - neigh_r%flux_bnd(1, 2))
            neigh_r%my(ng + 1) = neigh_r%my(ng + 1) + min_dt / neigh_r%dx * (current%flux_bnd(2, 3) - neigh_r%flux_bnd(1, 3))
            neigh_r%mz(ng + 1) = neigh_r%mz(ng + 1) + min_dt / neigh_r%dx * (current%flux_bnd(2, 4) - neigh_r%flux_bnd(1, 4))
            neigh_r%bx(ng + 1) = neigh_r%bx(ng + 1) + min_dt / neigh_r%dx * (current%flux_bnd(2, 5) - neigh_r%flux_bnd(1, 5))
            neigh_r%by(ng + 1) = neigh_r%by(ng + 1) + min_dt / neigh_r%dx * (current%flux_bnd(2, 6) - neigh_r%flux_bnd(1, 6))
            neigh_r%bz(ng + 1) = neigh_r%bz(ng + 1) + min_dt / neigh_r%dx * (current%flux_bnd(2, 7) - neigh_r%flux_bnd(1, 7))
            neigh_r%en(ng + 1) = neigh_r%en(ng + 1) + min_dt / neigh_r%dx * (current%flux_bnd(2, 8) - neigh_r%flux_bnd(1, 8))
         else
            write (*, *) 'Error from subroutine block_boundary:'
            write (*, *) ' Current level has no right neighbors and it is not a boundary block.'
         endif
      endif
   else
      call bnd_free_x (current%ro, ng, current%nx, 'r')
      call bnd_free_x (current%mx, ng, current%nx, 'r')
      call bnd_free_x (current%my, ng, current%nx, 'r')
      call bnd_free_x (current%mz, ng, current%nx, 'r')
      call bnd_free_x (current%bx, ng, current%nx, 'r')
      call bnd_free_x (current%by, ng, current%nx, 'r')
      call bnd_free_x (current%bz, ng, current%nx, 'r')
      call bnd_free_x (current%en, ng, current%nx, 'r')
   endif

   do i_block = 1, current%nc
      child => current%child(i_block)
      call block_boundary (child)
   enddo

   return
end subroutine
